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Abstract 

We address the nonadiabatic quantum dynamics of macrosystems with several coupled electronic 
states, taking into account the possibility of multi-state conical intersections. The general situation 
of an arbitrary number of states and arbitrary number of nuclear degrees of freedom (modes) is 
considered. The macrosystem is decomposed into a system part carrying a few, strongly coupled 
modes, and an environment, comprising the vast number of remaining modes. By successively 
transforming the modes of the environment, a hierarchy of effective Hamiltonians for the environ- 
ment is constructed. Each effective Hamiltonian depends on a reduced number of effective modes, 
which carry cumulative effects. By considering the system's Hamiltonian along with a few members 
of the hierarchy, it is shown mathematically by a moment analysis that the quantum dynamics 
of the entire macrosystem can be numerically exactly computed on a given time-scale. The time 
scale wanted defines the number of effective Hamiltonians to be included. The contribution of the 
environment to the quantum dynamics of the macrosystem translates into a sequential coupling 
of effective modes. The wavefunction of the macrosystem is known in the full space of modes, 
allowing for the evaluation of observables such as the time-dependent individual excitation along 
modes of interest, as well a spectra and electronic-population dynamics. 
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I. INTRODUCTION 



aradigm : 



Within the last twenty years, so-called conical intersections have slowly emerged as a 
aradigm for ultrafast nonadiabatic processes occurring in polyatomic systems 

Conical intersections are widely spread, actually omnipresent, types of intersect- 
ing potential-energy surfaces, leading to a complete breakdown of the Born-Oppenheimer 
approximation. The electronic and nuclear motions are strongly coupled in the vicinity of 
the intersection, which thus provides a very efficient pathway for the ultrafast decay of the 
excited electronic state (s), typically on the femtosecond time scale. 

Most of the studies reported to date exhibit intersections involving two electronic states. 
A notable exception is the case of triply-degenerate states due to symmetry requirement, 
as found, for instance, in the methyl cation [9]. More recent contributions also highlight 
accidental, i.e., occurring at nonsymmetric geometries, simultaneous degeneracy of three 



electronic states. The latter has been found, for instance, in the ethyl 



10(, allyl and 



pyrazolyl 121] radicals, in cytosine [13|] and in malonaldehyde Il5[ |. Within a given 
electronic-state manifold, cases of intersections between different pairs of states have also 
been reported, for instance, in C2H [igI and in the benzene radical cation jl7]. Of course, 
two-state and three-state intersections can be present within the same manifold of electronic 



states, see, e.g. 



1^. 



In this work, we aim at studying the multi-state non-adiabatic quantum dynamics in 
large polyatomic systems, i.e., those involving a large number of nuclear degrees of freedom 
(modes). We refer to them as " macrosystems" . These include, for instance, an impurity 
in a solid, a chromophore in a protein pocket, or a large isolated molecular system. The 
abundance of conical intersections grows with the dimensionality of the system and the 
density of electronic state in the energy domain considered. In large polyatomic systems, 
two-state intersections are known to be more common than avoided crossing situations 
Q. This togeth. With the fact that t.n.ee-state i.te.ect,o„. occur already ,n 

relatively small molecular species, let us suspect that multi-state intersections should be 
rather common and that several states are very likely to strongly interact nonadiabatically 
in macrosystems. We shall provide a theoretical framework to study the quantum dynamics 
in such situations. A special, but particularly important situation is given by the multi-mode 



Jahn- Teller effect 



19 



20|. 
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Often, the dynamics of macrosystems can be viewed to be dominated by a " system" part 
comprising a few strongly coupled modes only. Then, the large number of remaining modes 
is seen as an "environment". This environment may play an important role by modifying 
the quantum dynamics provided by the system alone. Indeed, the environmental impact 
on the system's dynamics has been highlighted in numerous situations, and is particularly 
important in conical intersection situations due to the high sensitivity of the latter to even 



small perturbations 
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271]. The impact of the environment should 



therefore be included in a realistic treatment of the macrosystem's dynamics. 

Nowadays, very powerfuU methods exist to treat numerically exactly the quantum dynam- 
ics of molecular systems of moderate size. We think^Jn particular, of the multiconfiguration 



time- dependent 



Hartree (MCTDH) method |28|, 



30 



3l| . and to its multi-layer extension 



(ML-MCTDH) 32|, |33| . MCTDH is able to treat the quantum dynamics of systems involving 



20-30 modes in two-state-intersection situations and of somewhat less modes if more states 
are involved. When even more modes are involved, as is typically the case in macrosys- 
tems, one natural approach is to use approximate quantum dynamical schemes which can 
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35( 1 . However, the full, 



account for conical intersection situations, see, for instance, Refs. 
numerically exact, quantum treatment of the dynamics is out of reach. 

Recently, another strategy was proposed. The aim is to construct reduced models which 
can account for, at least, the dominant effects of the environment upon the system. The 
limited number of modes in such models allows for a numerically exact treatment of the 
quantum dynamics. Thus, quantum aspects such as interferences, geometric phase effects 
361], etc, are naturally taken into account. Such an approach has been proposed for two- 
state conical intersections in macrosystems 25|. There, it has been shown that the use 
of three effective environmental modes only -together with the system's modes- suffice to 
calculate accurately the band shape and short-time dynamics of the entire macrosystem. 
Detailled analysis of this effective-mode theory along with numerical applications can be 



found in Refs. 
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38| . Precursors of this approach were derived more than twenty 



years ago for the Jahn- Teller effect 39 



40|, 



41 



42 



43, 



441 ]. This approach allows to split 



the environment into two parts: (i) a primary set of three effective modes which couples 
to the system's modes and carry the environmental effect on a short-time scale, and (ii) a 
"residual environment" which couples only to the effective modes and becomes important 
at later times. Importantly, the three effective modes are constructed from all the modes 
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of the environment and carry the cumulative effects of the latter on the short-time scale. If 
time scales beyond short-times are under interest, the use of these three effective modes is, 
however, not sufficient, and one has to take into account the residual environment. This can 
be done by constructing additional sets of effective modes, as recently proposed in Ref. {45]. 
Based on extensive numerical examples, it has been highlighted that the systematic use of 
additional effective modes allows to calculate accurately the quantum dynamics for longer 
and longer times. In this vein, a related extension of the effective-mode theory has been 



used to analyze exciton dissociation in semiconducting polymers [46l. l47l| . 

In this theoretical paper, we extend these recent findings to the general case of multi- 
state intersections in macrosystems, and analyze the dynamical properties of the proposed 
approach. We shall provide a detailled scheme for the construction of successive sets of 
effective modes in the situation where any number of electronic states are coupled. More 
specifically, we consider the general case of an arbitrary number of coupled electronic states 
and an arbitrary number of environmental modes. It will be shown mathematically that 
the use of the system's modes augmented by a limited number of sets of effective modes 
suffice to calculate, numerically exactly, the quantum dynamics of the entire macrosystem 
on a given time-scale. 

The paper is constructed as follows. In section II, we discuss the Hamiltonian used to 
describe the macrosystem with several coupled electronic states. In section III, we detail the 
construction of the sets of effective modes, which allows us to decompose the Hamiltonian 
of the environment into a hierarchy of effective Hamiltonians. In section IV, the dynamical 
properties of this hierarchy are analyzed and discussed. Section V concludes. 



II. THE HAMILTONIAN OF THE MACROSYSTEM 

We start with the Hamiltonian H of the macrosystem. H describes coupled electronic 
states in a diabatic representation. From now on, the "hat" symbol refers to N x N matrices 
in the electronic space. We define H as follows: 

H = Hs + Hb, (1) 

with a "system" Hamiltonian given by 

Hs = {hf,}, t,j = l,...,N, (2) 
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where the elements hf ^ are Hamiltonians depending on the Ns nuclear degrees of freedom 
(modes) of the system. These Hamiltonians /if^- are not further specified here, since all 
the developments which will follow do not depend on their particular form. Apart form 
containing the kinetic energy operators for the system's modes, they can, for instance, 
include the full (diabatic) potential energy surfaces obtained from ab-initio data for these 
modes. In general, Hs can couple all the N electronic states. In this work, the rest of the 
Hamiltonian, Hb-, which we shall call "Hamiltonian of the environment", is described by 
the linear vibronic coupling (LVC) model This model corresponds to the first term in 
a Taylor expansion of the actual potential energy surfaces at a given point in configuration 
space (reference geometry), in terms of so-called normal modes. This well established model 



has been successfully used to treat the quantum dynamics o 



molecular systems involving 
The Hamiltonian Hb is 



several coupled electronic states, see, for instance, Refs 
thus given by: 

Hb = hoi + {h^}, (3) 

with 

Nb 

ho = J2^{pl + 4), (4) 

k=l 

Nb 

= Y.'^k'^^k, for t,j = l,...,N (5) 

k=l 

where Nb is the number of modes of the environment and 1 the unit matrix in the electronic 
space. The Xk are position operators and the pk the corresponding momentum operators 
expressed using mass- and frequency-weighted quantities. The diagonal elements of Hb 
consist of harmonic oscillators, ho, augmented by state- and mode- dependent shifts given by 
/i^j. The off-diagonal element hfj couples the electronic states i and j. All the A''^ modes 
of the environment can, in principle, couple all the electronic states. 

Equations ([I]) to ([5]) determine the Hamiltonian of our macrosystem which contains 
Ns + Nb modes and involves A^ coupled electronic states. The above definition of the 
Hamiltonian of the macrosystem constitutes a natural extension of the 2-state Hamiltonian 



used as starting 



tems 
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26 



271, 



3oint in closely related former works on conical intersections in macrosys- 



37|. 



Let us discuss the Hamiltonian introduced above for A^ coupled electronic states. In 
Eq. we have isolated a "system" part from the rest ("environment") of the Hamiltonian. 



This separation is not always needed, but is usefull at least in two important situations, 
(i) When we aim at studying the impact of an environment on a "system", it is obviously 
convenient to use the proposed form of the total Hamiltonian. (ii) When some modes of 
the macrosystem (presumably a few, see below) cannot be reasonably treated with the LVC 
model, they can be treated as being part of a "system", which, in our approach, is not 
restricted to any kind of model. Of course, if all the modes of the macrosystem can be 
satisfactorily described by the LVC model, we can include all of them in Hb- In this case, 
Hs does not describe any dynamics and reduces to a diagonal matrix with the energy of 
the electronic states at the chosen reference geometry as elements. Note that the modes 
entering Hs and Hb do not couple directly, but do couple indirectly via the electronic 
subsystem. This is easily recognized by noticing that {hfj} in Eq. provides an explicit 
system-environment interaction. 

We make some further remarks on the LVC model. First, this model is valid for symmetry 
enforced as well as accidental intersections of potential energy surfaces. In this work, we 
even allow for the possibility that each environmental mode may contribute simultaneously 
to the diagonal as well as off-diagonal parts of the Hamiltonian. This allows one to account 
for cases where the environment lacks any symmetry. Of course, the LVC model can also 



describe all standard special situations, among which is the Jahn- Teller effect |19l. l20l|. The 
LVC model can be thought of as the minimal general model for the dynamics of intersecting 
potential energy surfaces. 

We suppose throughout this paper that the number of modes of the system is not too 
large, so that one can compute the quantum dynamics provided by Hs- It is the potentially 
very large number of environmental modes which renders the computation of the quantum 
dynamics of the entire macrosystem impossible. In the following, we shall introduce sys- 
tematic transformations of the Nb modes of the environment which will allow us to exactly 
decompose Hb into a hierarchy of effective Hamiltonians. These effective Hamiltonians, 
which will depend on a reduced number of effective modes, will be shown in Sec IV to 
translate into a sequential description of the dynamics of the entire macrosystem. Thus, 
depending on the time-scale under interest, only the first members of the hierarchy of effec- 
tive Hamiltonians will be needed to compute the dynamics of the entire macrosystem. Since 
each effective Hamiltonian will depend only on a reduced number of effective modes, the 
dynamics can be computed numerically (if, of course, the number of effective Hamiltonians 
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taken into account is not too large.) This approach has been recently derived for the 2-state 



case and shown by a numerical example to produce very accurate results for the quantum 
dynamics, see Ref. |45| (see also a related perspective in Refs. In the followine; 

sections we construct the hierarchy of effective Hamiltonians for coupled electronic states 
and study its dynamical properties. 



III. CONSTRUCTION OF THE HIERARCHY OF EFFECTIVE HAMILTONI- 
ANS 

We shall introduce successive orthonormal transformations of the modes of the environ- 
ment which will allow us to build the hierarchy of effective Hamiltonians. The system part 
of the total Hamiltonian will not be affected by these transformations of the environmental 
modes. The first step in the construction of the hierarchy will allow us to split Hb into two 
parts: 

HB = Hl+Hrl. (6) 

The first part, Hi, which constitutes the first member of the hierarchy of effective Hamil- 
tonians, will be shown to take account of all the couplings between the electronic states 
due to the environment and will thus play a key role. Importantly, Hi will depend only on 
a limited set of effective modes, constructed from all the modes of the environment. The 
second or remaining part of Hb, called Hri, will be diagonal in the electronic space. 

In a second step, we will iteratively construct additional sets of effective modes out of 
the remaining modes of Hri and obtain a hierarchy of effective Hamiltonians Hm'- 

Hrl = H2 + Hr2 
Hr2 = H3 + Hr3 

Hrm = + -f^rm+l; (7) 

where, in each step, the next member of the hierarchy, Hm+i, is constructed from the former 
remaining part Hrm- 
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A. Construction of the first eff"ective Hamiltonian 



We construct in the following the first member of the hierarchy of effective Hamiltonians, 
Hi. We shall give a detailed scheme for this construction; a similar scheme will be used 
in the subsequent construction of the rest of the hierarchy. The content of this subsection 
can be viewed as an extention of the results for the 2-state case to the general 

A^-state case of interest in this paper. To construct the first member of the hierarchy, we 
shall identify effective modes in Hb, and introduce accordingly an appropriate orthonormal 
transformation of all the Nb modes of the environment. 



1. Introduction of the effective modes 



Inspecting Hb, we see that we can define effective modes. Since Hb is hermitian, there 



exists a maximal number of A/" = A^(A^ + l)/2 independent elements hfj for i,j 
see Eqs. (E]) and We write them in terms of Af effective modes Xi, I = 1, . . 

hf^^ = R^'^'^Xi, l = {t-l)N + j-t + l, 

with effective coupling constants R^'^'^^ defined by 




(8) 



(9) 



We assume that the number of environmental modes, Nb, is larger than AT; if it is not the 
case our approach is not usefull as such. 

It is convenient to use a matrix- and vector-notation. We introduce the column 
vector of the Nb initial modes x = (xi, . . . , Xtv^)"^, the vector of the J\f new modes, 
Xj^ = {Xi, . . . , Xj\f)'^ , and a matrix Vj^xNb- The bold-faced quantities indicate vectors 
and matrices in the configuration space, and the subscripts indicate the size of the matrix 
or vector (no subscript means the full space, A''^). We define T^xAfs as: 



V, 



ATxNb 



(1,1) 
^(1.1) 

,(1.2) 

n 

5(1.2) 



(JV-1,JV) 

!H 

K(JV-l,iV) 
„{Af.^V) 

21 

k(^.JV) 



(1,1) 

''Nb 
5(1.1) 
(1.2) 
B 



\ 



5(1.2) 



(N-1,N) 
K(JV-l.iV) 

K(JV.iV) 



(10) 
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We straightforwardly have 

XaT = VxxNbX- (11) 

The X;, / = 1, . . . , A/", are already normalized, but are not orthogonal to each other. We 
orthogonalize them by using a standard method, e.g., the Gram-Schmidt orthogonalization 
procedure, and call the corresponding orthogonal modes Xi, I = l,...,J\f and the Af- 
dimensional vector they form = (Xi, . . . , X_\f)'^. We thus have 

X^ = LVxA/'^AT, (12) 

with Uj\fxj\f being a, M ^ M matrix which orthogonalizes the modes X;, / = 1, ...,A/'. 
Combining Eqs. (ITT]) and (I12p we arrive at 

= TxxNb^ (13) 

where 

TmxNb = Uj^xAfVAfxNB- (14) 

This defines a set of Af orthonormal effective modes, X;, / = !,... , A/", constructed from all 
the original modes x^. A; = 1, . . . , Nb- 

By applying the transformation Tj^^Nb to the original modes of the environment, the 
terms hfj of Hb can be expressed as 



/.5 = «:(^'^)VirPx, (15) 



with the coefficients K^^'^^^ given by 



1=1 



(ij) 

r\j 



k=l 



where the tik are the elements of TX/'xtyb- Note that the k[^'^'^ satisfy the normalization 
conditions 'Yl^=i{^i^'^''Y — 1 J — 1,---,X. The coefficients Ki^'^"* represent how the 
contributions to the effective coupling constant k'^*'-'^ are distributed among the M effective 
modes. 

In Eq. (IH]), each hf^ depends on a single (non-orthogonal) effective mode. After the 
orthonormalization of these modes, each element hf^ can depend on all the M orthonormal 
effective modes defined in Eq. (fT3|) . Importantly, by applying the transformation TVxats, we 



have transformed Eq. ([5]) into Eq. (fT5|) : the hfj depend now only on the effective modes, 
instead of the Nb original modes. It is already evident now that these effective modes will 
play a crucial role in the dynamics of the macrosystem, since they suffice to represent all 
the coupling terms hfj present in the Hamiltonian of the environment. 

We have constructed, so far. A/" effective modes which we call the primary effective modes. 
The full Hamiltonian of the environment contains Nb modes and, therefore, in order to fully 
describe Hb we need to identify Nb — M additional modes. This amounts to construct 
a full, A'^B- dimensional orthogonal transformation matrix. We call this matrix T^^\ The 
superscript (1) means that this is the first transformation related to the construction of Hi + 
Hri, see Eq. ([6]), and that subsequent transformations will be used later when constructing 
the other members of the hierarchy. Defining the complete set of new modes as X = 
(Xi, . . . , Xtv^)^, the matrix T^^^ is obviously such that 

X = T^^^x. (17) 

In fact, it is easily recognized that the matrix ^Vxtys contains the Af first rows of the com- 
plete orthonormal transformation matrix T^^\ These rows determine the primary effective 
modes. The remaining Nb — M rows of T^^^ can be chosen in many ways as long as T^^'^ 
is orthonormal (see Sec. III. A. 2 for a particularly appealing choice). Obviously, the trans- 
formation in Eq. f|T7|l leads to the elements hf^ given by Eq. f|T5l) . To fully transform the 
total Hamiltonian of the environment, we have to transform the diagonal part h^l as well, 
see Eq. ([3]). This gives 

Nb q Nb 

^0 = E y + + E ^«'(^;^^' + XiXv), (18) 
1=1 l<l'=l 

with Pi the momentum associated with Xi and 

Nb Nb 

= ''^^^ktfk, dill =''^^LJktlktl'k, (19) 

k=l k=l 

where the tik are the elements of full transformation matrix T^^\ Eqs. ( IT5|) and (|T8l) give 
Hb in the new, complete set of modes. 

We are now in the position to isolate the part of the transformed Hb which contains only 
the A/" primary effective modes and to obtain the desired form of Hb 

HB = Hi + Hri, (20) 
10 



where the Hamiltonian Hi containing only the M effective modes is the first member of the 
hierarchy and reads 

1=1 1=1 
N 

+ du'{PiPv + XiXi,)i. (21) 
i<i'=i 

The residual part Hyi takes on the diagonal form in the electronic space 

Nb ^ 

Hri = E Y^Pi+xfn 

l=N+l 

Nb 

+ E duiPiPv + XiX,)i 
1<V=M+1 
N Nb 

+ E E du>{PiPv + XiXi,)i. (22) 

1=1 i'=j\f+i 

It is worth noting that since the transformation is complete and orthonormal, the new 
form of Hb in Eq. (!20l) is completely equivalent to the original one given by Eq. ([3]): these 
two versions of the Hamiltonian describe obviously exactly the same physical problem. In 
the new version of the Hamiltonian, Hi depends on three terms. The first one consists of 
harmonic oscillators, the second one couples the electronic states, and the third one contains 
additional bilinear kinetic and potential terms which couple the M primary effective modes 
among themselves. The remaining part, Hri-, also contains three terms. The first (harmonic 
oscillators) and second (bilinear kinetic and potential coupling terms) ones depend only 
on the Nb — M remaining modes. Importantly, none of these modes participates in the 
coupling between the electronic states! The third term is made of bilinear kinetic and 
potential couplings between the M primary effective modes and the remaining modes. All 
terms of Hri are diagonal in the electronic space. Thus, the transformation T*^^) decomposes 
Hb into two parts, one. Hi, which does couple the electronic states and one, Hri-, which 
does not. Hi depends on M effective modes only. For the 2-state case, we recover that 



we need three effective modes 



25 



37(1 . For the 3-state case, 6 effective modes and for the 
4-state problem, 10 effective modes are needed, and so on. We remind that due to specific 
properties of some macrosystems, e.g., those of high symmetry, one may need less than M 
primary effective modes to construct Hi, but never more. The number of effective modes is 
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indeed equal to the number of linearly independent elements hfj in the original Hamiltonian 
of the environment Hb, see Eq. (|5]). Various special cases of the 2-state case have been 
discussed in Ref. S^], and can now be easily extended to the present context. Among 
others, a particularly important case is the multi-state multi-mode Jahn- Teller problem. 



2. A unique choice for the orientation of the effective modes 

To derive Hb = Hi + Hri it was unnecessary to explicitly specify two "quantities": 
(i) the particular choice of the orthogonalization matrix U^xjv, see Eq. f lT^ . and (ii) the 
choice of the Nb — M remaining rows of the full orthonormal transformation matrix T*^^) in 
Eq. (fT3l) . The final results for Hi and Hri in Eqs. (!2Tl) - (!22l) are valid whatever these choices 
for these quantities are. The freedom of choice is related to (i) the orientation of the M 
primary effective modes which span Hi on one hand, and (ii) the orientation of the Nb — M 
remaining modes on the other hand. 

All choices of orientation within the two subspaces of modes leads to mathematically 
equivalent results. However, some particular choices can lead to usefull simplifications of 
Hi and if^i- These simplifications concern the bilinear kinetic and potential coupling terms 
among the effective modes within Hi on the one hand, and among all the modes within H^i 
on the other hand. We want to stress that these simplifications are not needed to construct 
the remaining members of the hierarchy (see below), but give rise to working equations in a 
more closed form. Moreover, this form is particularly amenable to interpretation. 

We now introduce a particular orientation of the M modes entering Hi. Indeed, it 
is possible to exactly remove the third term of Eq. (pTI) . z.e., to remove all the bilinear 
coupling terms among the M primary effective modes. To proceed, we first notice that 
once the M modes have been constructed as detailled above, one can rotate them among 
themselves without changing the space they span, and thus without altering the impact 
of Hi. The particular orientation of interest is the one which diagonalizes the M ^ M 
matrix with elements Yl!k=i^ktiktvk-, 1,1' = I, . . . ,A/', see Eq. (fT9l) . The eigenvalues of this 
matrix are now the frequencies of the rotated modes, and the set of eigenvectors provides 
the A/" X A/" rotation matrix from the initial set of effective modes to the new, rotated one. 
The procedure is fully described in Appendix E of Ref. [s^] for the 2-state case, and is easily 
extended to the A^-state situation discussed in this paper. This procedure corresponds to a 
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unique choice of the orientation of the TV primary modes, or, equivalently, to a unique choice 
of the orthogonahzation matrix C/aTxat, see Eq. f fT2|) . 

The same procedure, but now for the remaining Nb — M modes, exactly removes the 
second term of H^i in Eq. fl22]) by appropriately choosing the orientation of these modes. 
This procedure fixes the remaining Nb — A/" rows of the transformation T^^^ in Eq. (1171) . As 
a result of this procedure the Hamiltonian takes on the simplified form 



Hb = Hi + Hrl, (23) 

Hi = J2 Y^P^ + + {^-^ Yl (24) 
1=1 1=1 

Nb ^ 

Hrl = E y(Pf + Xf)i 

l=Af+l 

M Nb 

+ duiPiPi' + XiX,% (25) 

1=1 i'=j\f+i 

where we have kept the same notation as in Eqs. (!20l) - (!22l) . although all quantities refer to 
the modes rotated as described above. Note that the A/" effective modes entering Hi are now 
only coupled through the electronic subsystem, and that in H^i the remaining Nb —M modes 
are decoupled from each other. The equations fl25]) - fl2S]) constitute our working equations. 

Apart form the fact that the mathematical form of Eqs. is simpler than the 

one of Eqs. (!20|) - (!22|) . which is of interest by itself and has numerical advantages, the use of 
the above unique choice of orientation of the M effective modes and of the residual modes 
is motivated by the following physical arguments. The first argument concerns Hi which is 
the only part of Hb which couples the electronic states directly. In Eq. fl2T|) . the primary 
effective modes are coupled among themselves in two ways, indirectly via the electronic 
subspace, and directly by the bilinear kinetic and potential coupling terms. By using the 
proposed orientation, the direct couplings among the primary effective modes are eliminated. 
As a particular consequence, the kinetic energy operator of the modes entering Hi takes on 
the usual form, and all the "unpleasant" momentum-space couplings disappear. This has 
an important implication. Now, if the off-diagonal elements of the matrix Hi go to zero, 
i.e., if there is no electronic coupling, the primary effective modes are no more coupled to 
each other as is the case in the usual picture of the LVC model. If, on the other hand, other 
choices of orientation are used, see also the related discussion for the 2-state situation in 
Refs. (271, l37| , the modes remain coupled to each other. 
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The second argument concerns the residual part Hrii and is related to the distributions 
of the couplings between the primary effective modes and the residual modes as a function 
of the frequencies of the latter. Here, again, the unique orientation of the A''^ — M residual 
modes chosen suppresses the bilinear couplings among them and leads to the usual form 
of the kinetic energy within the subspace of residual modes. This choice also leads to a 
unique set of M distributions of the bilinear couplings dni between the primary and residual 
effective modes (one distribution for each primary mode /, see Eq. These distributions 

explicitly reflect all the coupling between the primary and residual effective modes as the 
residual modes are not coupled among themselves. This provides a particular physical 
meaning to the above distributions, in addition to the simpler mathematical form of the 
Hamiltonian obtained by the unique choice of modes. 



B. Construction of higher members of the hierarchy 

In the preceding section we provided the first member of the hierarchy of effective Hamil- 
tonians. We shall now pursue the construction of the hierarchy. This is done by iteratively 
transforming the Nb — A/" modes of Hri which is the residual part of Hb, given by Eq. fl25l) . 
These successive transformations will allow us to decompose Hri into a hierarchy of effective 
Hamiltonians H2 + + . . . a,s follows: 

Hrm ~ Hjyi + H^^j^\. (26) 

The underlying idea is to construct additional sets of effective modes which successively carry 
cumulative effects of the respective residual part of the environment. The residual part Hrm 
of the environment shrinks as the hierarchy proceeds. This idea has been recently used in 



the construction of a hierarchy of effective Hamiltonians for the 2-state case in Ref. |45 |. 
We also refer to the recent work of Tamura et al. for a closely related construction 

of a hierarchy of 2-state effective Hamiltonians in the context of exciton dissociation in 
semiconducting polymers. The hierarchy of effective modes derived previously for the 2- 
state case is extended here to the A^-state case. 
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1. The second member of the hierarchy 



To construct H2 and Hr2 out of Hri-, we employ a very similar approach to the one used 
above to define Hi and Hri out of Hb- The approach consists of three steps: 

(i) identify additional collective (effective) modes in H^i, 

(ii) orthonormalize these effective modes, 

(iii) rotate these orthonormalized modes to obtain a more convenient mathematical and 
physical form of the equations. 

The first step is achieved by appropriately rewriting the term which couples Hri to Hi, i.e., 
the last term of Eq. (1251) . We keep for the moment only the position-dependent part and 
write 

Af Nb at rl 

J2 E du'X,X, = J2diXi £ ^X,, (27) 

1=1 l'=Af+l 1=1 l'=M+l ' 

where we define the effective coupling constants di by 

Nb 

d'i= / = 1,...,A^. (28) 

In Eq. (1271) we readily recognize M additional, normalized, effective modes defined by Xj^j^i = 
'^^Zj\f+iidu' / di)Xii , 1 = 1,... ,M . Equivalently, in a matrix notation: 

{Xjsf+i, X2mY' = WR-Ar)(^Ar+i, • • • , XNg__^f)'^ (29) 



with 



^ di,Af+i/di ■■■ di^2M/di ^ 



(30) 



ydj^^j^^i/d^ ■■■ dj^,2MldN j 
The number of these additional modes is of course equal to the number of primary effective 
modes. As seen in Eq. (|27j) . these modes are coupled to the primary modes. They are 
obviously orthonormal to the M primary effective modes, but are not orthogonal among 
themselves. We orthogonalize them using an identical procedure as the one given in section 
III.A.l, and call these M orthonormalized modes X^, A; = A/" + 1, . . . , iM: 

(A'aT+I, • • • , ^2nY = t^l/xAr(-^A^+l! • • • ! ^Nb-mY (31) 
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with Uj^l being the orthogonahzation matrix. Using Eqs. ([29]) and ([3l]) we obtain 

(X_v+i, . . . , X2mY' = ^Arx(Ars-Ar)(^Ar+l, • • • , Xj^g^j^)'^ (32) 

where 

T-'(2) _ T7-(2) 1/(2) foo\ 

The matrix corresponds to the first M rows of a full, (A^'^ — M) x (A^'^ — A/"), 

orthonormal transformation matrix T^^-*. As it was the case for T*^^) in the construction 
of the first effective Hamiltonian, the remaining rows of T(2) are chosen such as T^^^ is 
orthonormal. We finally obtain our vector of new orthonormalized effectives modes 

{Xj^+i, . . . , Xnb)^ = T^'^\Xj^+i, . . . , XnbV- (34) 

Transforming accordingly Hri, we obtain 

Hrl = H2 + Hr2 (35) 

where H2 is the second member of the hierarchy of effective Hamiltonians and reads (we 
drop the double tilde on the modes for simphcity) 

'2 



k=J\f+l 



JV 2N 

+ Yl Aik{PiPk + XiXk)i, (36) 

1=1 k=M+l 



and the residual part now takes on the appearance 



k'=2M+l 

2N Nb 

+ duk'{PkPk' + XuXk>)l. (37) 

k=N'+l fc'=2A/'+l 

Here, we have further used the unique choices of orientation of the A/" modes of H2, and of 
the corresponding residual Nb — 2A/' modes of Hr2 which simplify the equations as done in 
section III. A. 2 for Hi and H^i. The procedure to determine these orientations is exactly the 
same. The quantities Ai^, Vt^ and d^k' are easily determined by the analogous equations to 
those for the quantities for Hi, i.e., Eq. (fT6|) and (fT9l) . respectively. 
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We see from Eqs. fl5Bl) - (lTn) that only the M secondary effective modes entering H2 couple 
now to the M primary effective modes of Hi, whereas all the Nb — Af modes of Hri do so 
before the transformation. The Nb — 2A/' remaining modes of Hr2 couple to the secondary 
effective modes, and do not couple directly to the primary effective modes. Of course, since 
Hri is diagonal in the electronic space, H2 and Hr2 are also diagonal. 

2. The complete hierarchy 

Inspecting the new residual part Hr2 in Eq. fl37p . we immediately see that it is exactly 
of the same mathematical form as Hri in Eq- fl25|) . except that Hr2 depends only on the 
last Nb — 2Af modes which are the residual modes of Hi + H2- Consequently, we can use 
again our procedure described above to determine the third and higher members of the 
hierarchy of effective Hamiltonians. Eventually, this amounts to iteratively transforming 
the full Hamiltonian Hb of the environment to the sum 

HB = Hi + Y,Hm, (38) 

m=2 

where the primary effective Hamiltonian Hi is given by Eq. (12^ . and the higher members 
of the hierarchy Hm for m = 2, . . . , A'^v' all have the same formal structure and read 

mAf Q 

Hm= E -Y^P' + Xf')^ (39) 

l'={m-l)M+l 

(m-l)Ar mAf 

+ Yl AuiPiPv + XiX,)i. 

Except of Hi, all the H^ are diagonal in the electronic space. The number of sets of J\f 
effective modes is given by A'^v' = IntegerPart {Nb/M), with the integer part being here the 
smallest integer greater or equal to Nb/M. To simplify the following discussion, we shall 
call each set of N effective modes a multiplet of modes. Note that the last multiplet may 
include less modes than Af. 

Thus, Hb contains a total of A'^v' multiplets of modes. Hi, which couples the electronic 
states directly, depends on the first multiplet only. Each higher member of the hierarchy Hm 
couples to Hm-i and is diagonal in the electronic space. Consequently, while in the original 
Hb all the modes played formally the same role, after the successive transformations we get 
a sequential coupling of the multiplets of modes. 

17 



If one includes all the members of the hierarchy in a calculation on the macrosystem, of 
course, the exact results are recovered since all modes of the full Hamiltonian are taken into 
account. Generally, this is not possible in practice for a macrosystem because the number of 
environmental modes, A^^, or, equivalently, the number of multiplets Nj^, can be very large. 
In any practical calculation, we will truncate the Hamiltonian of the macrosystem. In this 
respect, the new form of the Hamiltonian Hb of the environment given by Eq. (l38!l . allows us 
to truncate the environment in a very systematic manner by resorting to the highest number 
of effective Hamiltonians one can or wants to afford for the problem at hand. It is indeed 
evident that truncating the hierarchy in this way is on a completely different level of quality 
than simply neglecting some of the environmental modes in the initial Hb- Each member 
of the hierarchy of effective Hamiltonians is constructed from all the modes of the residual 
environment, and carries cumulative effects represented by the effective coupling constants k 
and d. Furthermore, only the first multiplet of modes couples directly the electronic states, 
while in the original Hb all the modes of the environment couple these states. 

We shall reveal in the following the intimate relationship between the hierarchy of ef- 
fective Hamiltonians and the dynamical properties of the full macrosystem. The accuracy 
achieved when using a truncated hierarchy for computing the quantum dynamics of the 
entire macrosystem comprised of an arbitrary number of coupled electronic states and of an 
arbitrary number of modes is studied in the next section. 

IV. PROPERTIES OF THE HIERARCHY AND MOMENTS ANALYSIS 

We shall show in the following that a truncated hierarchy of effective Hamiltonians allows 
us to compute numerically exactly the quantum dynamics of the entire macrosystem on a 
given time-scale. The number of effective Hamiltonians included in the calculation will 
be shown to provide the corresponding time-scale. By including more and more effective 
Hamiltonians, the dynamics can be numerically exactly computed on a longer and longer 
time-scale. This central result is proven here by using an analysis of the moments of the full 
macrosystem. 
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A. Moments analysis 



We consider the autocorrelation function P{t) of the full macrosystem 



P(t) = (0|e-*^*|0) 



(40) 



where |0) is a column vector in the electronic space with components Tj|0), i = 1, N. |0) 
is the initial nuclear wavefunction in the initial - usually ground - electronic state and the 
are the transition matrix elements between this initial state and the manifold of coupled 
electronic states. Depending on the experiment performed, |rjp is the oscillator strength 
or ionization cross section for the state i. We suppose troughout this paper that the 
are (complex-valued) constants, i.e., are not position-dependent, and assume J2iLi = 
1. We suppose also that the initial nuclear wavefunction |0) can be factorize as |05)|0b), 
where \0s) is the initial wavefunction of the system and \0b) the noninteracting ground- 
state wavefunction of the environment. Thus, \0b) is separable with respect to the original 
environmental modes, and remains separable after our series of orthonormal transformations 
which lead to the final set of effective modes for the fully transformed Hb- Note that 
this is true because we use mass- and frequency-weighted coordinates and momenta. For 
convenience, we write \0b) as a direct product of nuclear wavefunctions for each multiplet 
of effective modes, rather than for the individual effective modes. Importantly, the initial 
wavefunction of the system, \0s), is not restricted to any particular form, as is the case for 
the system's Hamiltonian Hs- In particular, |0s) can include nonseparable contributions 
with respect to the system's modes. Our initial nuclear wavefunction is thus given by 



where the subscripts label the multiplets of environmental modes (and not the modes them- 
selves), with 



for all m = 1, . . . , A/". We assume |0) to be normalized. 

The autocorrelation function P(t) measures the overlap between the initial wavefunction 
and that at later times t |48|, |49|. The moments of the Hamiltonian are obtained by a Taylor 



|0) = |05)|0i)|02)---|0^J, 



(41) 




(42) 
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expansion in time of the autocorrelation function 

= E ^M„, (43) 

n=0 

where is the n-th order moment of the full macrosystem and is given by: 

M„ = (0|if"|0). (44) 

As the time t increases, more and more moments will contribute to the dynamical evo- 
lution of the macrosystem. The spectrum of the macrosystem corresponds to the Fourier 
transform of the autocorrelation function. Thus, the moments M„ are also connected to 
the spectral properties of the macrosystem. Mq is the total intensity. Mi gives the energy 
location, M2 is related to the width and M3 to the main asymmetry of the spectrum, etc. 
It has been shown for the 2-state case that Hs + Hi reproduces exactly the moments of the 



macrosystem up to and including M3 [25|, l37|] and this has opened the field of effective modes 
for multi-mode conical intersections. Recently, the hierarchy of effective Hamiltonians has 
been introduced for the 2-state case 
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46|. We extend this result to the full hierarchy with 
N electronic states and formulate a theorem: 

Consider a macrosystem with N coupled electronic states described by the Hamiltonian 
H = Hs + Hb of Eq. (Cp. The system's Hamiltonian augmented by the n first members of 
the hierarchy of effective Hamiltonians, Hs + Ylm=i^m, suffices to reproduce exactly all the 
moments with k <2n + 1 of the entire macrosystem. 

Below we present a rigorous proof of this theorem. Before doing so, we mention that 
while writing this paper we learned that Tamura et al. |47|] discussed the moments in the 
2-state situation without considering a system part Hs and came to the same conclusion for 
this special case. 

We start the proof with the fully transformed Hamiltonian, Eq. (!38|) . which contains Nj^ 
multiplets of effective modes build from all the Nb modes of the environment, plus the Ns 
modes of the system 

H = Hs+Y,Hm (45) 

ni=l 

and rearrange as 

H=Y,ifi'm + Smi). (46) 
m=l 



20 



where 

ni = Hs + Hi~Sil n^ = Hm-S^lWm>2. (47) 

Sm is the zero-point energy of the m-th multiplet of modes and we also introduce S = 
X]m=i ^m, which corresponds to the zero-point energy of the full Hamiltonian of the envi- 
ronment. Note that, for convenience, we include the system part in Tii in Eq. (147|) . i.e., only 
Til couples the electronic states. We stress that the Hamiltonian H in Eq. fH6|) is completely 
equivalent to the original Hamiltonian of the entire macrosystem given by Eq. ([1]), since the 
entire hierarchy as well as the system are considered. The moments M„ of the Hamiltonian 
H are thus exactly the moments of the entire macrosystem. From Eqs. ( Hl|) - fH^ we obtain 

N^f 

Mn = (0|i^"|0) = (0|(^7i:„ + ^i)"|0) (48) 

m=l 

Nat 



fc=0 ^ ^ m=l 

tQs-'M,, (50) 



where the coefficients (^) are the usual binomial coefficients, and the partial moments Aik 
are given by 

^fc = (0|(^7^™)l0). (51) 

m=l 

To evaluate the moment M„ we have to evaluate the partial moments Aik for all k < n. 

Before proceeding with the calculation of the moments, let us give two usefuU properties 
of the Hm- 

(PI) n^\0) = Vm > 2 (52) 
(P2) \nk,ni]=0 if l^k±l (53) 

The property (PI) is easily proven using Eqs. (147|) . (!39|) and (14T!) . It is also straightforward 
to prove (P2): Tik depends on the multiplets k and k — 1, and Tii on the multiplets / and 
/ — 1, see Eq. (15^ . We readily see that, if / > /c + 2 or if / < A; — 2 then and Hi do 
not depend on the same multiplets and thus commute. For = / it is trivial that they 
also commute. Recall that the Ti^ are diagonal in the electronic space for all A; > 1. The 
two properties (PI) and (P2) follow directly from the particular form of the initial nuclear 
wavefunction and of the Hamiltonians which compose the hierarchy. 
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We return to the evaluation of the moment Mn and evaluate the term M.k oi Eq. (150|) 
separately for odd and even powers of k which we denote for simplicity 21 and 2Z + 1, 
respectively. We write 

M21 = (0|(5^?i^)'(^7i:™)'|0) 

m=l m=l 

= m (54) 

Nat N^r 

m=l m=l 

= (55) 

where we define 

\L) = (5^7i:„)'|o) = Y.'^rnlLil). (56) 

m=l m=l 

Let us evaluate |1). Making use of the property (PI), one immediately finds 

|1) = J]7i™|0)=7i:i|0). (57) 



m=l 



The state 12) reads 



1 2) = J2nm\l) = Yl ^-^i|0) (58) 

m=l m=l 

= {A + n2)ni\o), (59) 

where the last equality makes use of (PI) and (P2): in the sum over m, according to (P2), 
all Tim with m > 3 commute with ?ii and can thus act directly on 1 0) , giving zero according 
to (PI). 

We now introduce 

Bi = Y,'Hm, (60) 

m=l 

where, of course, / is limited to Nj^f, the total number of effective Hamiltonians. The operator 
Bi contains all the Tim for m < /. Trivially, 

Bi+i = Bi + Hi+i, (61) 

and all -B; are hermitian since they are defined as sums of hermitian operators. We recall 
that the system's Hamiltonian is included in Tii, see Eq. (1471) . and thus included in all the 
operators Bi, I = 1, . . . , Nj^. With this new notation we can write Eqs. ( 1571) and ( l59l) as 

|1) = 4|Q); |2) = 52|1). (62) 
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We now prove by induction that \ l + 1 ) = A+ilO- Let us suppose by hypothesis that 

\l) = Bi\tll), (63) 

and evaluate \ l + 1 ). Using Eq. (156!) . one finds 

Nat 

11+1) = (64) 

m=l 

Nat 

= (5^7t:™)A|/^), (65) 

m=l 

where the sum runs over all the multiplets, i.e., the full Hamiltonian is considered. However, 
in the summation, the Ti^ with m > / + 2 do not contribute to | / + 1 ) because Bi contains 
only the Tik with k < I and because of the properties (PI) and (P2). Consequently, only 
the terms Hm with m < / + 1 in the sum contribute and, therefore, 

/+i 

= iJ2^rn)Bl\l_^) (66) 

m=l 
l+l 



:J2^mm (67) 



m=l 



= Bi+.IL), (68) 

where we have used Eq. f pUl) . Thus, if the hypothesis Eq. (jHS]) is true for |/), it is true for 
| / + 1 ). Finally, since we have already shown that the hypothesis is valid for Z = 1,2 (see 
Eq. (!62|) ). we have proven its validity for all /. Consequently, we immediately get: 

\L) = BiBi_,---B,B,\0). (69) 

With this result we can now evaluate A^2n and A42n+i from Eqs. (15^ and (1551) : 

M21 = {0\Bi--- Bi^iBfBi^i ■ ■ ■ fii |0) (70) 
M21+1 = {0\B, ■ --BiBi+^Bi ■ ■ -EilO). (71) 

To obtain our final result, one last step is required for M.2i+i- This quantity can be 
reduced to 

M 21+1 — 

= ai(A+H;+i)io = (iiAio (72) 
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because of 

ai^mlO = 0. (73) 

To prove the last equality, one needs to explicitly consider the state |/) in order to evaluate 
expectation values. Introducing explicitly the electronic basis, this state can be defined as 

|/) = (|/«),|/(^)),...,|/(^))f, (74) 

where the superscript (z) labels the component of the state |Z) in the corresponding electronic 
state i, and 

|/^'^) = l/S...)|Ow)---|0;v^), (75) 

where \fsi...i) depends on the modes of the system and on the first / multiplets of effective 
modes. Inspection of Eq. (!69|) shows that the system's modes and the multiplets k < I are 
inseparable in |/) as indicated by \fsi...i) while the multiplets k > I + 1 remain unaffected 
when constructing Hi+i depends only on the multiplets / + 1 and /, see Eqs. ( EUh and (H7I) . 
We write Hi+i as QiJf^i{Pl^-^^ + Xf^-^^) /2 — Sijf^i + di^i^i{PiPi^i + XiXi^i) where the indices refer 
to the A/" modes of the corresponding multiplets / and / + 1. Evaluating the expectation value 
(Ll'Hi+ill) using |/) in Eq. (17^ . we readily see that all terms vanish: the harmonic oscillator 
contribution is cancelled by i^j+i, and the bilinear part does not constribute because the 
expectation values of and P/+i in \Oi+i) are zero. This concludes the proof of Eq. (1721) . 
We now have 

M21 = (0|5i---5,_ifif5,_i---5i|0), (76) 
M21+1 = (0|fii ■ ■ ■ Bi_^BfBi_^ ■ ■ ■ B,\0) (77) 

which proves that only the system's Hamiltonian augmented by the members Tim with m < I 
of the hierarchy of effective Hamiltonians are needed to exactly evaluate Ai2i and M.21+1 
and all the lower-order partial moments M.k with k < 21. In turn, this implies that by 
using the truncated hierarchy Hs + Ylm=i we reproduce exactly all the moments Mk with 
A; < 2n + 1, see Eq. (!50|) . This proves the theorem and is the central result of this paper. 
As an immediate consequence one finds that each time one takes into account an additional 
member of the hierarchy, two additional moments of the macrosystem are recovered exactly. 
Further properties of the hierarchy are discuss below. 
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B. Discussion 



The hierarchy of effective Hamiltonians allows us to express the original Hamiltonian 
of the environment in terms of sequential couplings of multiplets of effective modes and is 
closely related to dynamical and spectral properties of the entire macrosystem. Using a 
given number of multiplets, the quantum dynamics on a given time-scale is numerically 
exactly reproduced, whatever the number of environmental modes is. This paves the way for 
studying truly large macrosystems involving a manifold of coupled electronic states. 

On the very-short time-scale, only the first member of the hierarchy plays a role. It 
the only member of the hierarchy which couples directly the electronic states as does the 
system's Hamiltonian Hs- All the impact of the environment onto the system on this very- 
short time-scale is contained in this first member which determines the three first moments of 
the macrosystem. Then, at later times also the second member becomes relevant. Its effect 
is to spread the vibrational energy whiting each electronic state. There is no direct energy 
transfer between the system and the effective modes of the second member of the hierarchy. 
This transfer is mediated by the first member. At even later times the third multiplet comes 
also into play. This multiplet now spreads further the energy of the macrosystem, and so on. 
Whenever a new multiplet comes into play, the energy is further spread into more directions 
in the environment. 



1. Autocorrelation function and band shape of macrosystems 

Due to the large number of environmental modes and the presence of several coupled 
electronic states, the autocorrelation function, which measures the overlap between the initial 
wavefunction and the time-evolving one, decreases usually very rapidly. This is a very well 
known feature for 2-state conical intersections, where the autocorrelation function typically 
decays in the 10-100 femtosecond time-scale, and subsequently exhibits some oscillations 
0, [J. These oscillations are usually of small amplitude because of the multimode nature of 
the dynamics which spread the wavefunction into many directions and lowers substantially 
the overlap with the initial one. The initial decay plays a crucial role. After this initial decay, 
the autocorrelation function typically possesses small values. This is of central interest here, 
since a few members of the hierarchy suffice to reproduce numerically exactly the short-time 
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dynamics, and thus the ultrafast initial decay of the autocorrelation function. 

By Fourier transforming the autocorrelation function, one obtains the spectrum of the 
macrosystem. The above remarks on the autocorrelation function can be translated in 
the frequency domain as follows. The fast initial decay of the autocorrelation function 
determines the band shape of the spectrum, which is thus accurately reproduced by using 
a limited number of multiplets. The autocorrelation function at later times after the initial 
decay determines the fine structures of the spectra, carving the band-shape of the spectra. 
The typically small value of the autocorrelation function after the initial decay translates in 
small carving of the band-shape. Again, this is well known for 2-state conical intersections, 
where the spectra are often so dense that individual spectral lines are difficult to resolve 
even for relatively small molecular species {3, Q] • By using the proposed approach, one can 
thus obtain accurate band shapes of truly large macrosystems by employing a few multiplets 
of modes only. These band shapes can be compared to experimental results. 

2. How many moments of the Hamiltonian are independent? 

We have proven above that the truncated hierarchy Hs + J2m=i '^rn reproduces exactly 
all the moments of the entire macrosystem up to M2„+i. Now, for the sake of discussion, 
let us assume that we consider the full hierarchy, i.e., n = Nj^. One then gets exactly all 
the moments up to M2Nj^+i- On the other hand, since the hierarchy is complete we have 
taken account of the full Hamiltonian, and this implies that all the moments are exact, and 
not only the 2Nj^ + 2 first ones (starting with Mq) . We may conclude that there are only 
2Nj\f + 2 "independent" moments of the full macrosystem. If one determines correctly the 
first 2Nj^ + 2 moments using our systematic approach, it follows immediately that all the 
higher moments are also correct. Let us consider an example. Assume there are Nb =24 
environmental modes in a =3-state problem. Each multiplet of modes contains A/" =6 
modes (sextet), and there are A^v" =4 sextets in total. This means that there are only 10 
"independent" moments for this problem: Mq to Mg. All the higher-order moments are 
automatically reproduced exactly if these first 10 moments are exact. We mention that the 
"independence" of the moments is to be understood with respect to the environment. The 
moments also depend on the system part, which is fully included in our theory, and for 
which all the moments can be "independent" depending on the complexity of Hs- 
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We stress that the property of the moments found above is a consequence of the form 
of the hierarchy of effective Hamihonians, and the definition of the multiplets of effective 
modes. If one uses instead the original form of the Hamiltonian of the environment, all the 
original modes must be included in order to compute even the first non-trivial moment M2. 
And, of course, once all the modes are included in the calculation, all the moments of the 
Hamiltonian are correctly reproduced. This is, however, of little help since we cannot include 
all these modes in any realistic calculation. Consequently, even the width of the spectrum 
(related to M2) cannot be reproduced exactly when using the original Hamiltonian of the 
macrosystem. 

3. Truncation of the hierarchy vs. truncation of the moment expansion of the autocorrelation 
function 

The following remark is of fundamental importance. Computing the quantum dynamics 
provided by the system's Hamiltonian augmented by the hierarchy truncated at the member 
n < Njij-, one reproduces exactly the 2n + 2 first moments (Mq to M2„+i) of the full Hamil- 
tonian of the macrosystem. We stress that the quantum dynamics provided by a truncated 
hierarchy does not lead whatsoever to a truncation of the Taylor expansion of the autocorre- 
lation function - which would have been dramatic for the dynamics on a longer time-scale. 
This can be easily seen. First, let us recall the definition of the moments, see Eq. (jS]), using 
the truncated hierarchy at the order n instead of the full Hamiltonian. The first moment 
which is not exactly reproduced by this approximate Hamiltonian is M2n+2 and reads 

n 

M2.+2 = + Yl ^™)'"^'|0)- (78) 

m=l 

Obviously, this moment does not vanish. Only the contributions arising from the first 
neglected term Hn+i of the hierarchy are missing in Eq. (1781) in order to make this moment 
exact too. Similar arguments hold for all other higher-order moments. Therefore, we not 
only obtain exactly the 2n + 2 first moments, but also important contributions to all the 
higher-order moments. The results on the dynamics at times longer than those reproduced 
accurately by the 2n + 2 first moments are approximate, but do contain the impact of 
important contributions of all the higher moments. 
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4- The wavefunction of the macrosystem 



As time proceeds, the propagation of the exact wavefunction depends on an increasingly 
growing number of members of the hierarchy of effective Hamiltonians. Let us call i„ the time 
up to which the dynamics is accurately reproduced by employing only the n first members of 
the hierarchy. In order to propagate the wavefunction, we consider in a dynamical calculation 
only the part of the wavefunction which depends on the first n multiplcts we want to account 
for. The remaining part needs not to be explicitly propagated and evolves trivially as 
free oscillators. Consequently, even when using a truncated hierarchy, we know the full- 
dimensional wavefunction of the entire macrosystem and all quantum observables can be, in 
principle, evaluated. 

We do not have a mathematical criterium to define a priori the value of t„ up to which 
the dynamics is accurately reproduced by employing the truncated hierarchy of effective 
Hamiltonians. Note that this value will depend also on the system's parameters. However, 
it is possible to estimate this value numerically by using the following procedure. We know 
that when we add a member of the hierarchy, two more moments of the full Hamiltonian are 
exactly reproduced. Thus, if we compare the numerical results (autocorrelation functions) 
computed using n members to those obtained by using n + 1 members, the results will start 
to deviate at t„. In fact, it is sufficient to add a single mode of the multiplet number n + 1, 
and compare the two autocorrelation functions obtained with and without this additional 
mode in order to get an estimate of tn- The wavefunction computed for times later than tn 
will be of diminishing quality as time increases, but will nevertheless account for important 
contributions of the exact wavefunction. This is easily seen from the discussion in point 3 
of this section. 

5. Effective modes vs. normal modes 

Our last point in this discussion concerns the link between the complete set of Nb effective 
modes and the original normal modes. To construct the effective modes, we have applied 
successive orthonormal transformations (rotations) of the original normal modes, see Sec. 
HI. These rotations can, of course, be concatenated into a single one, which we call T, and 
which is nothing else than the product of all the successive rotation matrices. Applied to the 
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original modes a;, this single rotation leads directly to the complete set of effective modes 
X which enter the hierarchy of effective Hamiltonians: X — Tx. It is obvious that, by 
using the transposed transformation T'^, we recover the original modes out of the effective 
modes: x = T^X. 

Since the full-dimensional wavefunction is known numerically exactly up to t„, see the 
last point discussed above, highly accurate expectation values can be computed in the full 
space of effective modes. On the other hand, knowing the transformation matrix from 
this space to the original space of normal modes, one can compute expectation values of 
quantities defined by the latter modes. Examples are the time-dependent energy gained by 
some normal modes of intersct and the average position of a particular normal mode. Such 
quantities are of help to interpret experimental results, where discussions in terms of normal 
modes are widely used. 

V. CONCLUSION 

The Hamiltonian of the macrosystem is decomposed into a system part and an environ- 
ment part. The latter is described by the LVC model and the former is comprised of a 
few modes which are not restricted to any model and usually include the modes which can- 
not be satisfactorily described by the LVC model. It is the huge number of environmental 
modes which exclude a direct, numerically exact calculation of the quantum dynamics of 
the macrosystem. 

Starting from the original Hamiltonian of the environment we derived a hierarchy of 
effective Hamiltonians. Each member of the hierarchy of effective Hamiltonians has been 
shown to depend on a single multiplet of effective modes. The number of effective modes 
within each multiplet depends on the number of coupled electronic states. The first member 
of the hierarchy plays a key role since it contains all the direct coupling among these elec- 
tronic states due to the environment. Only this first member participates directly, together 
with the system part, in the coupling of the electronic states. The higher members of the 
hierarchy are all diagonal in the electronic space, i.e., they do not couple the manifold of 
electronic states. 

Furthermore, except of the first member, each member of the hierarchy is coupled only 
to the "former" member by bilinear kinetic and potential terms. As a consequence, the 
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original Hamiltonian of the environment where all the modes formally play a similar role, is 
replaced by a hierarchy of Hamiltonians which exhibits a sequential coupling of multiplets of 
modes. This sequential coupling immediately translates to a hierarchical description of the 
quantum dynamics of the full macrosystem. We have proven, by analyzing the moments of 
the autocorrelation function of the entire macrosystem, that each member of the hierarchy 
comes into play at a different time. The truncated hierarchy at the order n suffices to 
reproduce exactly the first 2n + 2 moments of the entire macrosystem. The higher moments 
are also reproduced but only approximately. In this way one obtains numerically exactly the 
dynamics of the macrosystem up to a finite time tn which grows as n increases. The total 
wavefunction of the macrosystem is known in the full space of vibrational modes. Thus, 
apart from the autocorrelation functions, spectra, and the time-evolving populations of the 
electronic states, we can also evaluate other quantities related to the wavepacket, like, for 
instance, the individual excitation along all normal modes. 

This work extends the 2-state effective Hamiltonian formulation for short-time dynamics 



251] to the more general case of coupled electronic states and lifts the restriction to 
short-time dynamics. We expect this work to be valuable in studying the highly complex 
behavior of electronically excited macrosystems involving multi-state intersections and a 
large number of nuclear modes. The use of a truncated hierarchy, i.e., of a reduced number 
of effective modes for the environment, can allow for a numerically exact treatment of the 
quantum dynamics in such macrosystems. Finally, we mention that the hierarchy of effective 
Hamiltonians is certainly of interest also for the use of approximate dynamical methods. 
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